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ABSTRACT 

We study the spatial fluctuations in the hydrogen ionizing background in the epoch follow- 
ing reionization (z ~ 5-6). The rapid decrease with redshift in the photon mean free path 
(m.f.p.), combined with the clustering of increasingly rare ionizing sources, can result in a 
very inhomogenous ionizing background during this epoch. We systematically investigate the 
probability density functions (PDFs) and power spectra of ionizing flux, by varying several pa- 
rameters such as the m.f.p., minimum halo mass capable of hosting stars, and halo duty cycle. 
In order to be versatile, we make use of analytic, semi-numeric and numeric approaches. Our 
models show that the ionizing background indeed has sizable fluctuations during this epoch, 
with the PDFs being a factor of few wide at half of the maximum likelihood. The clustering of 
sources dominates the width of the PDFs, so analytic models must take large-scale clustering 
into account. The distributions also show marked asymmetries, with a high-value tail set by 
clustering on small scales, and a shorter low-value tail which is set by the mean free path. 
The power spectrum of the ionizing background is much more sensitive to source properties 
than the PDF and can be well-understood analytically with a framework similar to the halo 
model (usually used to describe dark matter clustering). Nevertheless, we find that Lya forest 
spectra are extremely insensitive to the details of the UVB, despite marked differences in the 
PDFs and power spectra of our various ionizing backgrounds. Assuming a uniform ionizing 
background only underestimates the value of the mean ionization rate, (r^), inferred from 
the Lya forest by a few percent. Instead, analysis of the Lya forest is dominated by the un- 
certainties in the density field. Thus, our results justify the common assumption of a uniform 
ionizing background in Lya forest analysis even during this epoch. 
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1 INTRODUCTION 

The metagalactic hydrogen-ionizing ultraviolet background (UVB) 
is probably the most influential, pervasive radiation field through- 
out most of cosmic history. Much investigation has gone into under- 
standing its properties and its impact on the intergalactic medium 
(IGM) (see Mei ksinl |2007| for a recent review). At low redshifts 
(z < 3^1), the ionizing photon mean free path (m.f.p.) is very 
large, and so each gas element effectively "sees" many sources of 
ionizing radiation. A uniform UVB seems to be an accurate ap- 
proximation in this regime, with spatial fluctuations contributing 
only a few percent of the mean dZuoll 19921: iFardal & Shulll l 1993) : 
Gnedin & Hamilton! 120021 : iMeiksin & White! 120041 : ICrofd 12004 ; 



Bolton et alj2006l) . 



However, the m.f.p. decreases rapidly with redshift, with ex- 
trapolations predicting that it should appro ach values of A m f p ~ 
30-40 comoving Mpc by z ~ 5-6 (e.g. IStorrie-Lombardi et al.l 
1 19941: lMiralda-Escudej|2003l) . Such modest values of the m.f.p. at 
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z > 5 can induce significant fluctuations in the UVB, sourced by 
cosmic variance in the number of ionizing sources inside each A m f p 
horizon. Additionally, the 1/r 2 flux attenuation, combined with the 
clustering of sources, can int roduce add i tional scatter in the UVB 
on scales smaller than A mfp feud 1 19921 : [Meiksin & Whitell2004h . 
These fluctuations h ave been largely ignored i n high-red shift analy- 
ses (t hough see §5 in lBolton & Haehnelt 2007 and §4 in lFeng et al.l 

book) . 

Here we undertake a systematic study of the post-reionization 
metagalactic ionizing background at z ~ 5-6. We vary source 
prescriptions and A m f p in an attempt to infer general conclusions 
about the theoretical properties and observational importance of 
spatial fluctuations in the UVB. We focus on the simplest statis- 
tics, the probability density function (PDF) and power spectrum, 
which are most likely to be observed. Given the dearth of knowl- 
edge concerning the properties of absorbers and sources at high- z, 
we make use of a range of methods that include analytic, semi- 
numeric and numeric approaches in order to be comprehensive. We 
compare simple analytic models of these quantities to the numeri- 
cal results in order to shed light on the origin and importance of the 
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variations; we will show that analytic models are accurate in some, 
but not all, cases. 

This paper is organized as follows. In Sj2] we outline the vari- 
ous modeling tools we use to study the fluctuations in the ionizing 
background. In Sj3]and ij4]we present the PDFs and power spectra 
of various UVBs, respectively. In [|5] we study the impact of a spa- 
tially varying UVB on the Lya forest. Finally, in Sj6] we summarize 
our findings and offer conclusions. 

We quote all quantities in comoving units, with the exception 
of flux, and we denote proper units with a prefix 'p'. We adopt 
the background cosmological parameters (JIa, SIm, Qb, n, as, Ho) 
= (0.72, 0.28, 0.046, 0.96, 0.82, 70 km s~ * Mpc" 1 ), matching 
the five-year results of the WMAP satellite (Dunkle v" et alj|2008l ; 
iKomatsu et"ail2008h . 



2 GENERATING INHOMOGENEOUS FLUX FIELDS 

2.1 Semi-Numerical Simulations with Flux Tool 

We genera te our ionizing flux fields following the procedure de- 
scribed in Mesinger & Dijkstra (2008). We briefly outline the pro- 
cedure below. 

We begin with halo fields at z = 5 and z = 6 gen- 
erated with the semi-numerical simulation DexlvQ, which has 
been shown to reproduce the correct number density and cluster- 
ing properti es of halos well into the q uasi-linear and non-linear 
regim es dMesinger & Furlanetto! 120071 : Fig. 1 in IDiikstra et al.l 
120081 : Mesinger et al. in preparation). Our simulation boxes are 
150Mpc on a side, with a 150/1800 Mpc halo grid cell size. The 
velocity fields used to perturb the halo field were generated on a 
lower resolution 900 grid; thus our final halo field resolution is 
Ax = 150/900 = 0.17 Mpc. 

For each halo field, we create a corresponding UV flux field 
on a 150 3 grid (spatial resolution of 1 Mpc), by performing a halo 
mass/r 2 weighted sum. Specifically, we compute the flux of ioniz- 
ing photons (in units of ionizing photons s _1 pcm~ 2 ) with 

/(x,*) = £±i£ ]T Xi(DC) ^_ Xei °; e-'— 

Mi >M min ' 1 

(1) 

where x is the location of the cell of interest, Mi is the total halo 
mass, Xi is its location, and the factor of (1 + z) 2 converts the fac- 
tor |x — Xi| 2 from comoving into proper units. We also include a 
duty cycl^fl parameter through the random variable _X"i(DC), which 
has a value of 1 with likelihood DC or with likelihood 1 — DC. 
Finally, the factor eion in eq. Q} denotes the rate at which ioniz- 
ing photons are released into the IGM by a dark matter halo per 
unit mass. Unless stated otherwise, we assume a fiducial value of 
e ion (z) = 3.8 x 10 58 /DC [^^] ^ff ., which provides 
a good fit to the observed luminosity functions of Lya emitting 
galaxies (LAEs) dShimasaku etal] 120061; iKashikawaet al J 120061 ; 
IDiikstra et ai1l2007l ; IStark et al.ll2007l ; iMcOuinn et alj|20071) , and 

1 http://www.astro.princeton.edu/~mesinger/DexM/ 

2 Note that there are several definitions of the term "duty cycle" in the 
literature, some of which include this same parameter in the star-formation 
timescale of each halo (an assumption only valid in extremely idealized 
scenarios). We adopt the simple definition above of the duty cycle, namely 
the fraction of dark matter halos which are " on" (i.e. hosting active galaxies 
emitting ionizing photons into the IGM; c.f. lLee et alj2003) . 



the z = 6 Lyman Break galaxies (LBGs) l lBouwens et alj l2006). 
However, as we are mostly concerned with the shape of the flux 
PDF, where applicable we present results in units of the mean flux 
or intensity, J/ (J). 

In constructing the UVB, we also assume a minimum halo 
mass able to host stars, M m j n , and explore mainly two different 
values for this parameter: 1.6 x 10 8 and 1.4 x 10 9 Mq. The for- 
mer value corresponds to a virial temperature of 10 4 K at these red- 
shifts, and represents the regime of ine ffective feedback on atomi- 
cally cooled halos du ring reionization dMesinger & Diikstrai l2008l : 
lOkamoto et alj2008h . The later value represents the regime of very 
inefficient star formation inside such small hal os, perhaps due to 
strong radiative and/or mechanical feedback (e.g.lYepes et alj 19971 : 
IScannapieco et alj|2006l ; IPawlik & Schavel l2008). The exact value 
of the higher M m i n was also motivated by the fact that the number 
density of halos with mass greater than 1.4x10 Mq is roughly ten 
times smaller than the number density of halos with masses greater 
than 1.6 x 10 8 Mq at z ~ 5, thus allowing us to compare results 
at fixed source number density with values of DC = 1.0 and 0.1 in 
the two models, respectively. 

2.2 Hydrodynamical Simulations with Flux Tool 

Additionally, we generate a flux field corresponding to the z — 
5.71 source field f rom the co smological hydrodynamical simu- 
lation presented in iTrac et af. I d2008h (their z = 6 reionization 
model). This simulation is fixed-grid, 143 Mpc on a side, and in- 
cludes prescriptions for modelin g dark matte r, b aryons and ioniz - 
ing photons (for details see lTrac & Cen 2007| and lTrac et alj2008h . 
The density field was calculated on grid of 0.19 Mpc cells, which 
resolves the Jeans length in the mean density, ionized intergalactic 
medium (IGM) by a factor of few, and then smoothed to a cell size 
of 0.74 Mpc. Each halo's instantaneous star-formation rate (SFR) 
is proportional to the instantaneous gas accretion rate, and enters 
our eq. (TJ in place of the halo's mass, Mi, with the normalization 
adjusted accordingly so as to match the mean value of the UVB. 
We make use of this flux field, combined with the z — 5.71 gas 
density fields from the hydro-simulation, to generate more accurate 
mean Lya forest flux decrement statistics in Sj5] At these redshifts 
and scales, our semi-numerically generated density fields some- 
what over-predict the rare voids which dominate the this statistic 
(Mesinger et al,, in preparation). Furthermore, the ray-tracing algo- 
rithm from the numerical simulation over-predicts the fluctuations 
in the flux field, due to an insufficient number of ray splittings (Trac 
2008, private communication). This problem only becomes severe 
following reionization, but since this is the epoch we study here, 
we use this hybrid prescription in some of the results below (i.e. we 
do not use the radiative transfer field from the simulation). 

In Fig. Q] we present a 0.74 Mpc thick slice through the den- 
sity a nd flux fi e lds ba sed on this z = 5.71 simulation output 
from pTrac et al. I d2008l) . The flux fields were calculated assuming 
A m fp = 30 Mpc, and DC = 1. Even for such a moderately high 
choice of A m f p , we can qualitatively see that there is significant in- 
homogeneity in the flux fields. Furthermore, it is evident that the 
flux and density fields are highly correlated, an issue we will return 
to inSg] 

2.3 Analytic Model 

We complement these numerical techniques with a relatively sim- 
ple analytic model. While this approach cannot incorporate all of 
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Figure 1. Slices through the z = 5.71 gas density (left) and UVB (right) fields using the hydrodynamical simulation of lTrac et alj J2OO8) and our simple flux 
generation tool. Each slice is 143 Mpc on a side and 0.74 Mpc thick. The UV flux is generated using the star-formation rates obtained from the simulation and 
the prescription outlined in [ff] assuming A m f p = 30 Mpc, DC=1, and Af m ; n ~ 3 X 10 s Mq. 



the physics provided by the semi-numeric and numeric models, we 
will find it to be helpful in elucidating the physics of the back- 
ground radiation field. As in the remainder of this paper, we will 
focus on computing two statistical descriptions of the flux field: the 
PDF and the power spectrum. 



2.3.1 Flux PDF 

The PDF of the radiation background J, normalized to its mean 
value (J), can be computed exactly for randomly distributed 
sources, provided that we assume a constant attenuation length 
A m fp for ionizing photons much sm aller than the Hubble length 
(Zuo 1992k iMeiksin & Whit3l2004l) . In this case, the crucial pa- 
rameter is the number of sources within each attenuation length, 
No = (47r/3)rii An lfp , where m is the source number density, and 



dp{> X\N ) 
dX 



sNq I dx x(j>(x) Im G(sx) 



ds exp 



sX + sNq / dx xcj>(x)Re G(sx) 



where X = J/ J*, 4>(x) is the source luminosity function, in units 
of the mean luminosity and normalized to have unit integral, 

J* = W(47rA mfp ) 2 , 



G(t) 



dur (w)e 1 



(3) 



and u = e _T /r 2 . In this limit, the mean background is (X) = 
SN , or (J) oc No J* oc L*r . 

At the end of reionization, A m f p > 20 Mpc and (typically) 
No > 10 5 . Thus, the Poisson fluctuations in the galaxy counts are 
very small and one might expect the ionizing background to be ex- 
tremely uniform. However, equation ([2} assumes that the source 
counts within each attenuation volume are Poisson distributed 



around the fixed value No; in actuality, deterministic halo clus- 
tering causes fluctuations in the expected number of halos within 
each attenuation volume, with a fractional amplitude ~ 6cr(A m f p ) 
where b is the mean bias of the star-forming galaxies. For our fidu- 
cial model at z = 6, we find b ~ 2.5, which yields fractional 
variations of order ~ 2 for A m f p = 30 Mpc. 

To crudely account for this deterministic clustering, we ap- 
proximate the underlying density field as linear on scales ~ A m f p 
(so that the PDF of density, p(S), is a gaussian with variance 
o" 2 (A m f p ) and vary the expected number of sources within each 
attenuation volume according to linear theory, so that No — * 
No{5,X m i p ), which depends on the conditional mass function 
for a region with a fractional overdensity S on a scale A m f p 
l lLacev&Colell 19931) . Thus 



dp{> X) 
dX 



J d6p(5) 



dp(X\N [l + bS]) 
dX 



(4) 



We shall see that this is only a crude approximation, for two rea- 
sons. First, our prescription for the expected source density is 

(2) 

not terribly accurate, since the typical fractional fluctuation in the 
galaxy counts is not too much smaller than unity, and because we 
effectively assign the modified source density to all neighboring 
regions as well. More importantly, nonlinear clustering is actually 
quite important to the PDF. 



2.3.2 Flux Power Spectrum 

Similar techniques also allow one to compute the correlation func- 
tion or power sp ectrum of the intensity field due to randomly dis- 
tributed sources feuo & Phinnev|[T993h . However, in this case it is 
possible to include (linear) clustering exactly by using arguments 
similar to the "halo" model that describes dark matter and galaxy 
clustering. 
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The intensity field around a single galaxy of mass m is 

j( r \M) = ^nMiexp(-r/\ mip ) 



(5) 



(4 7rr )2 

In our calculations, the normalization factor is set as in the semi- 
numeric model. The total intensity field is simply the sum of the 
profiles around each galaxy in the Universe. By analogy to the halo 
model, the power spectrum can then be written as the sum of two 
terms, one arising from correlations within a single source's en- 
velope, and one arising from correlations betw e en the locations 
of tw o sources dScherrer & Be rtschinger Il99ll : ICoorav & Shetr] 
2002). These two terms both depend on the normalized Fourier 
transform of the single source profile, equation {5}: 

arctan(fcA m f p ) 



uj{k) = 



(6) 



fcA m f p 

which is independent of galaxy mass (because it only describes the 
shape, not the amplitude, of the intensity profile). Note that uj — > 1 
for fcA mfp < 1 and uj — > 7r/ (2fcA m f p ) for fcA m f p > 1. 

The "single-source" term, describing correlations within ran- 
domly distributed sources' profiles, is 

Pis{k) = \uj{k)\ 2 ^ dmn(m)( — ) , (7) 

where n(m) is the mass function of the sources (or the halo mass 
function in our calculations) and p ga i = n t m ga i is the mean mass 
density of galaxies (and m ga i is the mean mass of the ionizing 
sources). Note that the single-source correlation f unction can also 
be de rived fro m the two-point p robability density dZuo & Phinnevl 
1993; see also Furlanetto 2009). Furthermore, the amplitude is de- 
termined by the luminosity function, but the shape is entirely deter- 
mined by our imposed attenuation profile. 

Similarly, the term describing the correlations between two 
sources may be written 



P2s(k) = \uj(k)\ 2 b 2 P lin (k, 



dmn(m) ( — — 



(8) 

where Pn n (k, z) is the linear matter power spectrum. Again, this 
two-point function has a shape fixed by the source profiles, but an 
amplitude that depends on the nature of the sources. In fact, the 
term in square brackets is simply equal to unity; we have left it in 
only for comparison to the halo model, where u(k) is a function 
of m and so must be left inside the integral. Because \u.j\ 2 — » 
1 for small k, the two-source term traces the linear matter power 
spectrum on these scales, before being truncated at fcA m f p ~ 1. 

If we take a simple model in which all sources have a fixed 
mass, we find that 



E 2 ! 

Pis 



b 2 niPn n (k). 



(9) 



Thus the two-source term becomes more and more important as 
the number density of sources, their bias, and the underlying linear 
clustering increase: this is simply because all of these increase the 
importance of deterministic clustering relative to the correlations 
inherent to the flux profile itself. 



3 THE FLUX PDF 
3.1 Attenuation 

As discussed above, we account for flux attenuation through the 
factor e _r,/Amf p/(47rr 2 ) in eq. (TJ. Thus flux attenuation is taken 



into account in a spatially homogeneous way through the free pa- 
rameter A m f p . Although such an approach misses the obvious cor- 
relation of absorbers (LLSs and DLAs) with sources, it had sev- 
eral benefits over using the typical approximate radiative transfer 
algorithms in cosmological hydrodynamical simulations. Firstly, 
our approach is versatile, easy to calibrate and to run quickly for 
large-scale boxes. This allows us to easily explore the effects of 
varying A m f p , which is importan t, given the large uncertain ties in 
this value at these redshifts (e.g. Faucher- Giguer e et al.l2008l) . Sec- 
ondly, our formalism allows for easy comparisons with analytic 
models, and provides direct insight into the features of the UVB. Fi- 
nally, and most importantly, large-scale cosmological hydrodynam- 
ical simulations do not directly resolve LLSs and DLAs; thus they 
too must employ analytic prescriptions to populate their box with 
absorbers. We also not e that the mean free paths at these redshifts 
is so large (> 20 Mpc; | Storrie-Lombardi et al][T9 94; Mada u et al.l 
ll999l ; lMiralda-Escude1l2003l : lFaucher-Giguere et al]2008l) that flux 
received by each cell is sourced by tens of thousands of sources, 
likely averaging out detailed radiative transfer effects. We thus do 
not expect our results to be significantly affected by the fact that 
we neglect the details of radiative transfer and the clustering of ab- 
sorbers. 

Nevertheless, stochastic fluctuations in the absorber popula- 
tion will help to increase the variance in the radiation field: if LLSs 
are primarily responsible for absorption, each attenuation length 
will contain only a single absorber on average and the absorption 
length will vary quite strongly. However, in practice > 50% of the 
absorption is typically due to lower-column density absorbers (see, 
e.g., the Appendix to Furlanetto & Oh 200^), so this is probably 
not a large factor for us. 

In Figure [2] we plot the flux PDFs for a L — 100 Mpc box 
with 0.5 Mpc cells at z = 10, A mfp = 10 Mpc, and M min = 10 8 
Mq . We show the effects of truncating the sum over contributing 
sources in eq. (Q} at distances of lA m f p , 2A m f p , and 3A m f p (left to 
right solid curves). The distributions appear to converge at 2A m f p 
— * 3A m f p . The high end tails of the distribution are dominated by 
nearby ionizing sources and are fairly insensitive to the effective 
horizon (i.e. the high-value tails are roughly constant as more dis- 
tant sources are included in the UVB calculation). When plotted in 
unnormalized units (i.e. J instead of J/ (J)), this overlap is more 
pronounced and the leftmost (blue) curve shifts to even smaller val- 
ues. In contrast, the low-end tails of the PDFs correspond to loca- 
tions far away from ionizing sources and so are very sensitive to 
the inclusion of more distant sources, which roughly translates to 
changing the amplitude an effective homogeneous UVB. If sources 
were homogeneously distributed, truncating the flux distribution at 
distances greater than nA m f p would underestimate the ionizing flux 
by a factor of (1 — e _n ). The clustering of sources results in the 
more complicated behavior seen in Fig. [2] 

The dotted curve in Fig. [2] neglects the exponential attenua- 
tion term in eq. Q}, and instead only includes the contributions of 
sources at distances r < A m f p . This is analogous to assuming a 
step-function optical depth so that r = for distances r < A m f p 
and t = co at distances r > A m f p . Such a prescript ion can sim- 
plify computation in a nalytic models (e.g ., IZuol ll992f) and numer- 
ical calculations (e.g., iBolton et al. 2006) and seems to result in 
fairly accurate flux distributions, though it slightly overestimates 
the number of faint regions. This overestimate worsens somewhat 
if one integrates source s only out to lA m fp but also includes expo- 
nential attenuation (e.g. lSantos et alj2008t) . 

In Figure[3] we explore the impact of varying A m f p . Flux PDFs 
at z — 5 (solid curves) and z = 6 (dashed curves) are shown as- 
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Figure 2. Flux PDFs at z = 10, assuming A m f p = 10 Mpc, M m ; n = 
10 8 Mq and DC=1.0, calculated according to eq. (TJ but truncating the sum 
over contributing sources at distances of 1 A m f p , 2A m f p , and 3A m f p (left to 
right solid curves). The dotted curve assumes a step-function attenuation, 
essentially r = for distances r < A m f p and r = oo at distances r > 
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Figure 3. Flux PDFs at z = 5 (solid curves) and z = 6 (dashed curves) 
generated assuming M m ; n = 1.6 X 10 8 Mq and DC=0. 1. Starting from 
the left side of the figure, curves correspond to values of A m f p = 20, 30, 40 
Mpc, respectively. For comparison, the dot-dashed curve shows the squared 
overdensity, A 2 = (p/p) 2 , taken from the hydrodynamical simulation at 
z = 5.71 with cell size 0.74 Mpc. 



Figure 4. Flux PDFs at z = 10, assuming A m f p = 10 Mpc, M m i n = 
10 s Mq and DC=1.0. The dotted curve was generated assuming a step- 
function attenuation, essentially r = for distances r < A m f p and r = oo 
at distances r > A m f p . The dashed curve is obtained in the same manner, 
but with randomized source locations (i.e. ignoring clustering). The solid 
and dot-dashed curves correspond to the analytical model described in the 
text, including and ignoring clustering on scales > A m f p as described in 
equation @, respectively. 

sumingM min = 1.6x 10 8 M© andDC=0.1. Note that the evolution 
of structure from z = 10 in Fig.|2]to z — 5 and 6 in Fig.[3]preferen- 
tially suppresses the low-end tail of the PDFs as regions containing 
few ionizing sources become extremely rare. This results in nar- 
rower and more asymmetric distributions. 

As expected, increasing the m.f.p. also results in narrower 
flux distributions, as the gas cells see increasingly representative 
volumes of space inside their horizon. We note however that the 
^mf P = 40 Mpc distributions are still fairly broad. Understand- 
ably, when plotted in unnormalized units (i.e. J instead of J/ (J)), 
the high-value tails overlap as they are dominated by nearby halos. 
Perhaps surprisingly, we see very little variation of the flux PDF 
with redshift, over the range z — 5 to 6. 



3.2 Source Clustering 

The clustering of sources broadens t he flux PDF and creates a 
high-value tail in th e distribution (e.g. Mesi nger & Diikstrall20oj 
iDiikstra et al.l2008l) . How strong is this effect and how much do the 
properties of the underlying sources affect these trends? 

In Fig. [4] we plot flux PDFs at z = 10, assuming A m f p = 
10 Mpc. The dotted curve in the figure was generated from the 
semi-numeric calculation assuming a step-function attenuation at 
A m f p = 10 Mpc, shown to be an adequate approximation in Fig. 
[2] The dashed curve is obtained in the same manner, but with ran- 
domized source locations. We see that ignoring source clustering 
can severely underestimate the widths of the flux PDF. Even the 
shape of the curves is quite different. 

Fig. E] also compares the PDF to our analytic model. The dot- 
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dashed li ne use s the analog of equation but for step- wise attenu- 
ation (see IZuoll 19921 and lFurlanettdl2008n . again assuming Poisson 
distributed sources. It closely matches the semi-numeric results for 
the same assumptions. 

The solid line includes linear clustering over the scale A m f p as 
described in equation (0. The clustering length of M ~ Af m i n ~ 
10 8 Mq halos is ~ 2 Mpc at this redshift, so the clustering on scales 
larger than the m.f.p. should be well into the linear regime and accu- 
rately predicted by the analytic model. Obviously, this simple pre- 
scription provides a relatively poor match to the full semi-numeric 
results: although the large-scale clustering does broaden the distri- 
bution by a comparable amount, the shapes still disagree. 

Part of this difference is easy to understand and is a conse- 
quence of the analytic model's simple prescription for clustering: 
the model not only fixes the local attenuation volume to have a 
given density, but it also effectively assumes that surrounding re- 
gions have the same number of sources. This artificially suppresses 
the contribution from distant galaxies (which provide 2/3 of the 
mean ionizing background in these models). In practice, points in- 
side low density regions should be exposed to higher flux values 
from neighboring higher density regions. Such long-wavelength ef- 
fects are missed in the simple analytic treatment, resulting in the 
spuriously long small-value tail in the PDF. 

However, by the same logic our prescription places high- 
density regions near high-density neighbors, so one might expect 
us to overestimate the extent of the variation at high J as well - 
but in fact the prescription underestimates the width there as well. 
This is because, in contrast to the low- J points, the radiation field 
at such high levels is not dominated by the background of distant 
sources but by nearby clusters of objects, for which smaller-scale 
clustering and nonlinear effects are important. 

This is not particularly surprising in light of the analogy we 
have drawn to the halo model: although that simple prescription 
does a remarkably good job of matching the dark matter power 
spectrum, it cannot easily be used to predict the actual densi ty dis- 
tribution of the IGM. Indeed. IScherrer & Bertschi nger ( 199 1]) have 
shown that the PDF of the fractional deviation of J from the mean 
value, 1 + 8j = J/ (J), is 
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where £j (ri,r2, •••,!>) is thep-point correlation function of the 
radiation field. Thus the PDF depends on all higher order correla- 
tions at zero-lag, which are evidently strongly affected by the non- 
linear clustering on small scales. The simple addition of large-scale 
linear clustering is not sufficient to match the semi-numeric results. 

To see the impact of the source properties on the clustering 
signature in the flux PDF, in Fig. [5] we plot several flux PDFs as- 
suming Amfp = 30 Mpc. The dotted, short-dashed and solid curves 
were all generated at z = 5, but assuming various values of DC 
and M m ixx- These two parameters drive the correlation function of 
sources and are crucial in c onstraining the U V luminosity to halo 
mass relation at high-z (e.g. lLee et alj 2008). From Fig. [5] we see 
that the low-end cut-off is generally unaffected by these parame- 
ters. This is to be expected since, as we mention above, the low- 
end tail of the flux distribution corresponds to voids distant from 
UV sources, and are sensitive to the effective horizon A m f p (which 
is identical in all these curves). Thus, most of the clustering imprint 
is in the high-value tails of the PDFs. 

The figure also confirms that the widest flux PDF results from 
the model in which sources are the most biased and rare: having 
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Figure 5. Flux PDFs assuming A m f p = 30 Mpc. Dotted, short-dashed and 
solid curves correspond to z = 5 and (DC, M min /M ) = (1.0, 1.4 X 10 9 ), 
(0.1, 1.4 X 10 9 ), and (0.1, 1.6 X 10 s ), respectively. The long-dashed curve is 
generate d from the z = 5. 71 star-formation rate output from the numerical 
results of Trac et al. ( 2()08j), as described in the text. 

the highest Af m i n and lowest DC. As the duty cycle is increased, 
the high-value tail retracts. A similar trend is seen as one lowers 

M min . 

As mentioned above, the choice of Af m i n = 1.4 x 1O 9 M0 
was motivated by the fact that (DC, M min /M©) = (1.0, 1.4 x 10 9 ) 
and (0.1, 1.6 x 10 8 ) result in the same number density of sources, 
and thus are useful in isolating the imprint of each these two pa- 
rameters. Even in our simple model where the duty cycle is not 
a function of mass, these corresponding curves are quite similar. 
Thus it is unlikely that the flux PDF, even if it could be measured, 
could provide a useful constraint on these parameters. 

The long-dashed curve in Fig. [5] shows the flux PDF obtained 
from the z = 5.71 sta r-formation rate o utput from the hydro- 
dynamical simulation of iTrac et al.1 ((2008), as described in £|2,2| 
Despite differences in the star-formation prescription, this curve 
agrees well with the above trends. Note that the slight offset at low 
J is attributable to the fact that these fluxes were calculated includ- 
ing sources within 1.5A m f p , instead of within 2A m f p , as was the 
case for the semi-numerical boxes. 



4 THE FLUX POWER SPECTRUM 
4.1 Analytic Results 

We begin our discussion of the flux power spectrum by exam- 
ining the "halo model" analytic results (which we expect to be 
more successful than the PDF). Figure [6] shows some examples 
at z = 10; here A 2 (k,z) = k 3 /(2n 2 )P(k, z). The solid curve 
takes M min = 10 8 A/ Q , DC= 1, and A mtp = 10 Mpc; the dot- 
ted curves show the corresponding one- and two-source contribu- 
tions. Clearly these two terms are cleanly separated in scale for this 
set of parameters. The two-source term dominates on large scales; 
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Figure 6. Flux power spectra at z = 10 computed with our analytic 
model. The solid black line takes Af m in = 1O 8 M0, DC= 1, and A m f p = 
10 Mpc; the dotted curves show the corresponding one- and two-source 
contributions. The curve with solid squares shows the corresponding one- 
source power spectrum computed from the semi-numeric simulations. The 
long-dashed line takes DC= 0.1 instead, and the short-dashed line takes 
A mfp = 40 Mpc. The dot-dashed line takes M min = 1O 1O A40. 

for feA m fp "C 1, it traces the linear power spectrum with an am- 
plitude determined by the mean galaxy bias. It is then damped at 
fcAmfp ~ 1: the short-dashed curve, which takes A m f p = 40 Mpc 
shows how the turnover scale shifts with the attenuation length. 

On small scales, the one-source term dominates: it has A 2 S oc 
k 3 uj(k) 2 oc fc/A mfp . Thus its amplitude decreases rapidly with 
the attenuation length - although this dependence is hard to sep- 
arate from the dependence on overall source density. The dashed 
curve with solid squares shows the one-source term computed from 
the semi-numeric simulations (with sources randomly distributed, 
so that the two-source term disappears). It provides an excellent 
match to the analytic model, except on scales near k ~ 10 Mpc -1 , 
where the discrete grid cells affect the semi-numeric power spec- 
trum. 

The dot-dashed curve takes a larger A/ m i n . Note that the power 
spectrum increases with the halo bias on very large scales. Also, 
the decreased source number density strongly increases the impor- 
tance of the one-source term, giving a much less clean separation 
between the one-source and two-source regime. 

4.2 Numerical Results 

We now compare these analytic predictions to the semi-numerical 
simulations. For the numerical calculations, we use the conven- 
tion A 2 (k,z) =k 3 /(2n 2 V) (|<5./(k,z)| 2 ) fc . InFig.|7j we plot the 
flux power spectra at z — 6 for both the analytic (solid curves) 
and semi-numeric (dashed curves) models, generated assuming 
Mmin = 1.6 x 10 s M Q and DC=0.1. Curves correspond to val- 
ues of A m fp = 20, 30, 40 Mpc, from top to bottom. 

The agreement between the two calculations is quite impres- 



Figure 7. Flux power spectra at z = 6 from the analytic model (solid 
curves) and the semi-numeric simulations (dashed curves, with Poissonian 
error bars). All curves assume A/ m i n = 1.65* 10 s Mq and DC=0.; the three 
sets correspond to values of A m f P = 20, 30, 40 Mpc, from top to bottom. 
The dotted curve shows the semi-numeric result at z = 5 for A m f p = 
20 Mpc. 

sive, especially at large scales. The analytic model clearly shows 
the separation between the two-source and single-source terms; 
with the semi-numeric results the difference is less clear because 
the boxes are not large enough to sample the regime fcA m f p 2> 1, 
where A 2 traces the matter power spectrum. However, the ampli- 
tudes of the analytic and semi-numeric calculations match reason- 
ably well on these scales: the differences in the two-source regime 
are attributable to slightly different mean halo biases in the analytic 
and semi-numeric calculations. 

However, the semi-numeric results have a steeper slope for 
the one-source component; this is most likely due to nonlinear 
clustering, which will steepen the correlation function on small 
scales, since Fig. [6] shows that the two different methods agree 
very well for the one-source term when clustering is ignored. Inter- 
estingly, however, the crossover scale between the one-source and 
two-source regimes is accurately predicted by the analytic model. 

The three sets of curves in Fig. [7J show that the m.f.p. has lit- 
tle impact in the shape of the spectrum over the scales probed by 
the semi-numeric simulations: it primarily serves to decrease power 
fairly evenly on all scales with fcA m f p < 1. As described above, for 
smaller wave-numbers the attenuation length has no effect: effec- 
tively, a larger attenuation length causes more damping relative to 
the linear power spectrum. 

Finally, the dashed curve shows the flux power spectrum from 
the semi-numeric simulations at z = 5 and A m f p = 20 Mpc. By 
comparing Fig. [7] and Fig. [3] we see that the flux power spectra 
are more sensitive to redshift than the flux PDFs. The differences 
are similar for other values of A m f p , and in the analytic model: in 
all cases, the amplitude of the one-source component is somewhat 
larger at the lower redshift. 

This is easily understood from the analytic formalism. The 
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Figure 8. Flux power spectra analogous to the models shown in Fig. [3] 

two source term is proportional to b 2 Pn n (k, z); each of these terms 
evolves with redshi ft, but in t he product the evolution nearly can- 
cels at high redshift ( Oh 2001). The one source term is proportional 
to the second moment of the mass function (when integrated over 
mass - eq. [7jl, which increases toward lower redshifts as more and 
more massive halos are assembled. 

A similar argument holds for the power spectra presented in 
Fig. [8] where we vary the DC and M m in parameters. The curves 
correspond to the same models presented in Fig. [5] Here we see 
drastic differences on moderate to small scales, because models in 
which the sources are rare (i.e. with low DC and high M mln ) tran- 
sition to the single-source regime on larger scales. If observed, this 
statistic could be a powerful tool in constraining these parameters. 
The flux field has the benefit of encoding information about ioniz- 
ing sources which are too faint to be detected by galaxy surveys. 
Unfortunately, such flux power spectra will be nearly impossible 
to extract from the Lya forest, even given a very large number of 
sightlines (see Fig.QT]and associated discussion). 



5 LYMAN a FOREST 

The most obvious application of our flux field analysis is in the in- 
terpretation of the Lya forrest spectra following reionization. Thus 
far, almost all estimates of the photoionization rate per hydrogen 
atom, Ti2 = r/10~ 12 s _1 , obtained fr om the Lya for est have as- 
sumed a uniform UV background (e.g. lFan et al. 2006). However, 
in !j3]we showed that UV background actually has sizable fluctua- 
tions in this regime. These fluctuations could potentially introduce 
bias and stocasticity in the determinations of T12 from Lya forest 

atte nuation measuremen ts. 

Bolton & Haehneljl d2007l) briefly addressed this issue by in- 
cluding a similar UV background model to our semi-numeric 
approach. They confirm t he results of lower redshift studies 
dGnedin & Hamilton! [20021 : iMeiksin & White] 120041 : iBolton et al.l 
120061) that assuming a homogeneous UVB underestimates the mean 



Figure 9. Top panel: Mean ionization rate corresponding to an effective 
Lya (black curves) and Ly/3 (red c urves) GP optica l depth, for the z = 5.71 
hydrodynamical simulation box o f lTracetalJ fcoOS), assuming A m f p = 30 
Mpc. Dotted curves use the density field from the simulation, but assume a 
uniform UVB. Dashed curves assume a fluctuating UVB from our simula- 



tions, but convolved with the Miralda-Escude et al. (20t)3) density distribu- 
tion. Dot-dashed curves use the fluctuating UVB and density field from the 
simulations, but randomly correlate the two. Finally, solid curves use the 
correlated fluctuating UVB field and underlying density field from the sim- 
ulations. Bottom panel: Ratio of the mean ionization rate for each model 
i and the "true" model shown as the solid curves in the top panel. The 
vertical dashed line demarcates the average value of r® ff ~ 3.6, roughly 
correspondin g to the valu e measured from the SDSS QSOs at this redshift 
(taken from lFan et al . 2006). 

ionization rate obtained from Lya forest absorption, (£12)- They 
find that this effect is of order ~1% - 10%, depending on A m f p 
and M m in (see their Fig. 3). However, for this analysis they used a 
single realization of a 30 h^ 1 Mpc simulation box. This box size 
is comparable to a single ion i zing photon m.f.p at these redshifts 
JStorrie-Lombardi et a l. 1994; Madau et alii 19991 : iMiralda-Escudd 
120031: Faucher-Gigue re et al. I l2008h and is several times smaller 
than the length of the GP troughs used to estimate (r^). Thus, 
their analysis was unable to fully take into account cosmic vari- 
ance. 

Using our large simulation boxes, we wish to systematically 
investigate this effect. As a first step, we do not create detailed 
mock absorption spectra. Instead, for a given flux field, we merely 
calculate the Lya optical depth, r a , for each cell, assuming ion- 
ization equilibrium. Then we average the flux decrements over the 
entire box to obtain a corresponding effective GP optical depth: 
Ta a = -ln(exp[— T a ]) x . 

In the top panel of Fig.|9]we plot (ri2) for each correspond- 
ing effective Lya (black curves) and Ly/3 (red curves) GP optical 
depth, for the z = 5.71 hydrodynamical simulation box, assum- 
ing A m f p = 30 Mpc. The cell size for this simulation is 0.74 Mpc. 
As mentioned above, the density field was calculated on a much 
finer grid of 0.19 Mpc, which resolves the Jeans length in the mean 
density, ionized IGM by a factor of few. The smoothed cell-size of 
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0.74 Mpc corresponds to the Jeans mass in the mean density ion- 
ized IGM, and is roughly comparable to the Doppler width of the 
Lya line at T = 10 4 K, as well as to the Keck ESI resolution at this 
redshift. Thus we expect our r cS <-» (ri2) mapping to be adequate 
[i ndeed, our results are cons istent with the mock spectra analysis 
of lBolton & Haehneltl fe007h l. but we stress that we did not create 
detailed mock spectra. 

The curves in Fig.[9]show the effects of various simplifications 
on the value of (T^) inferred from the r cff measured in the spectra. 
The most accurate model, using both our fluctuating flux fields and 
underlying density fields, is shown by the solid curves. By compar- 
ing the solid and dotted curves (or their ratios in the bottom panel), 
we see that assuming a uniform UVB makes very little impact on 
the estimated mean ionization rate per hydrogen atom. At low val- 
ues of r cff , the two curves are practically indistinguishable, with 
the difference increasing with r eff . At r° ff ~ 3. 6 (roughly corre - 
sponding to the observed value at this redshift, e.g. iFan et alj2 006). 
assuming a uniform UVB results in a ~ 5% underestimate in the 
value of (ri2). 

It is fairly easy to interpret this behavior by looking back at 
Fig-H] The optical depth approximately scales as r oc A 2 /T, where 
A = p/p is the gas density in units of the mean. Thus the opti- 
cal depth distribution can be approximately obtained by convolv- 
ing these two quantities. From Fig. [3] we see that the distribution 
of A 2 is much wider than all of the flux distributions. Thus a uni- 
form UVB, which treats the flux distributions as a delta function, 
is an accurate assumption. Because of the exponential attenuation 
involved, the r cfI statistic is sensitive to the low-end tail of the 
density distribution. As the ionizing background gets fainter (i.e. 
increasing r e ), the r e statistic is sensitive to an increasingly nar- 
row range in the low-end tail of the optical depth distribution. It is 
in this regime that the non-zero width of the ionizing background 
PDF becomes evident, as the convolution with A 2 is mostly evi- 
dent in smearing-out this low-end tail. Thus assuming a uniform 
UVB underestimates the true (F^), but only slightly, with this un- 
derestimate increasing with increasing r cft . 

Figure|9]also shows the impact of other common assumptions 
in the r eB *-> {F 12 } mapping. The dot-dashed curve is generated 
by using the same UVB and density fields as in the solid curve, but 
randomly pairing the two. This also results in a slight underesti- 
mate of {Fi2}, though a little more notable than the uniform UVB 
assumption. Since underdense pixels are preferentially farther away 
from most sources, where the UVB is lower (c.f. Fig. [TJ, not cor- 
relating the two fields results in underdense pixels receiving fluxes 
which are too strong, thus lowering (Ti2). Nevertheless, the effect 
is still small because of the relative width of the PDFs argument in 
the above paragraph. 

Finally, the most drastic effect, resulting in an underestimate 
of a factor of ~2 in our models, comes from using the analytic den- 
sity PDF from Miralda-Escude et al. ( 2000). This offset is normally 
attributed to errors in the fluctuating Gunn-Peterson approximation, 
which ignores effects such as the line wings, peculiar velocities, 
and blending. However, our numeric results make a similar approx- 
imation, so that cannot be the source of the error here. Instead, it 
may be differences between the assumed density distribution (gen- 
erated to fit numerical simulations at redshifts z ~ 2-4, using a 
simulation box of only 10 Mpc) and that of the simulation. At 
higher redshifts (i.e. r eff values), the transmission becomes domi- 
nated by rare voids JOh & Furlanettoll2005l : llidz et al.|[20o3) . and 
so cosmic variance is very important. Proper analysis requires large 
simulation boxes with accurate hydrodynamical treatment and di- 
rect higher-redshift outputs. However, if the semi-analytic method 
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Figure 10. r <-> (Tia) mappings using our semi-numerical simulations, 
assuming various flux models. In the top panel, solid curves correspond to 
z = 5, while dashed curves correspond to z = 6. 

is renormalized using a simulation (as is typical), the shapes of the 
curves are similar enough that the errors will be relatively small, so 
the approximation is still useful for understanding the evolution of 
(ri2). Furthermore, it is not clear how this offset is affected by the 
assumed thermal history of the gas, which sets the effective local 
smoothing scale. A detailed exploration of this is beyond the scope 
of this work, and merits its own study. 

In Fig. [10] we show the r off «-> (ri2) mappings using our 
semi-numeric simulations, assuming various flux models. Note that 
the semi-numerical density fields slightly overpredict the presence 
of very rare voids which dominate this statistic (Mesinger et al., in 
preparation), thereby over-predicting the transmission for a given 
UVB. It is remarkable how insensitive this mapping is to the details 
of the UVB, despite the notable differences in the flux PDFs and 
power spectra shown in the previous section. Again, this is due to 
the fact that the fluctuations in the UVB, and the differences among 
various models of the UVB, are overshadowed by the fluctuations 
in the density field. The only notable differences among the curves 
is due to the density field evolution from z — 6 to z — 5. 

We note that our r esults are in general agreement with 
iBolton & Haehneltl j2007l) . who found that assuming a homoge- 
neous UVB underestimates (Fi2), by a few percent over this range 
of A m f p values. However, they find a somewhat stronger depen- 
dence on A m f p then we do. This is most likely a result of cosmic 
variance attributable to the small simulation box they used to study 
this effect. 

Finally, we wish to quantify how easy is it to distinguish vari- 
ous UVB models from observations of the Lya forest. One can see 
from Fig.[8]that the shape of the flux power spectrum encodes infor- 
mation about the clustering properties of the sources. As mentioned 
above, the flux field (as potentially extracted from the Lya forest), 
has a great benefit in being able to probe the impact of ionizing 
sources below the detection thresholds of galaxy surveys. Hence, 
in Fig. [TTJwe show the optical depth power spectra at z = 5 and 
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Figure 11. Optical depth power spectra at z = 5 and A m f p = 30 Mpc, 
within the range 0. 1< r Q <6, approximately corresponding to the range 
probed by current instruments. 

A m f p = 30 Mpc for some of the same models as in Fig. [8] Unfor- 
tunately, the optical depth power spectra of these models are prac- 
tically indistinguishable from each other, despite the large differ- 
ences shown in Fig. [8] This result could have been foreshadowed 
from Fig. [TO] and is again due to the fact that the density field and 
its fluctuations swamp the UVB signal. 

However, we note that the density fluctuations decrease as red- 
shift increases, while fluctuations in the UVB increase. The UVB 
becomes more inhomogeneous not only through the evolution of 
A m f p and the bias of sources, but it is also modulated by the ion- 
ized bubbles throughout reionization (whenever that epoch may 
be). Thus inhomogeneities in the flux field become increasingly im- 
portant at higher redshifts, although whether this can be observed in 
the highly saturated forest at z > 6, where the density fluctuations 
themselves cause substantial fluctuations in r eff dLidz et alf2 006). 
is questionable. 



6 CONCLUSIONS 

We systematically study the expected fluctuations in the hydrogen 
ionizing background in the epoch following reionization (z ~ 5-6). 
Unlike at lower redshifts (z < 4), the UVB can be quite inhomoge- 
nous in this regime. This is due to the smaller mean free path of 
ionizing photons, as well as the clustering of increasingly rare ion- 
izing sources. 

We confirm that there is a sizable spread (the widths of the 
PDFs at half of the maximum likelihoods spanning factors of ~ 
2-4) and asymmetry in the distribution of ionizing fluxes in a cos- 
mological volume. Expected values of the m.f.p. are large enough 
(Amfp > 20 Mpc), and sources are ubiquitous enough, that the PDF 
of the ionizing background is asymmetric, with a tail extending 
to higher values. This high-value tail is set by clustering on small 
scales, and is insensitive to the contributions of distant sources (i.e. 



Amfp) ( Me singer & Diikstra|[2008h . We also find that the PDFs are 
most sensitive to the m.f.p., with higher values of A m f p increas- 
ingly truncating the l ow-value tails. Furthermo re, we note that the 
analytic formalism of iMeiksin & White! (2004) severely underesti- 
mates the width of the flux PDF at this epoch, since it ignores the 
clustering of sources. 

The power spectra of the ionizing background are even more 
sensitive to the details of the source model. We have shown that an 
approach analogous to the halo model, in which contributions from 
the proximity zones of single sources (which dominate fluctuations 
on small scales and are independent of the large-scale galaxy distri- 
bution) and multiple sources (which dominate on large scales, and 
which depend on the large-scale clustering of the ionizing sources) 
are separated, provides a good match to the simulated results, al- 
though nonlinear clustering provides even more small-scale power 
than the analytic model predicts. The single-source component is 
particularly sensitive to the number density of galaxies; unfortu- 
nately, these power spectra will be nearly impossible to observe 
with the Lya forest. 

We also model the imprint of various UVBs in the Lya forest 
at z ~ 5-6. We find that the Lya forest spectra are extremely insen- 
sitive to the details of the UVB, despite marked differences in the 
PDFs and power spectra of our various UVBs. This is attributable 
to the fact that the Lya optical depth scales as r tx A 2 /T, and 
the squared overdensity distribution, /(A 2 ), is much wider than 
all of our r distributions. In fact, even the extreme scenario of a 
uniform background (i.e. a delta function F distribution), only un- 
derestimates the value of (ri2) inferred from the Lya forrest by a 
few percent. Thus our results justify the common assumption of a 
uniform UVB in Lya forest analysis. 

We also confirm that accurate modeling of the density field is 
very important in the r cff «-» (r^) relation. We find significant 
differences (fact or of ~ 2) in this relation among the numerical 
and extrapolated Miralda-Escude et al. (2000) density distributions 
(though we stress again that we did not create detailed mock spec- 
tra in this work). At higher redshifts, the r cff statistic increasingly 
depends on a narrower range in the low-value density PDF. Large- 
scale hydrodynamical modeling of the statistics of these rare voids 
is therefore e ssen tial in interpreting the high-redshift Lya forest. 
IPawlik et all J2009h have recently studied the IGM density field 
during reionization, and do find significant differences from the 
Miral da-Escude et al.l (2000) fitting formula, although their simu- 
lation boxes are much too small to study the rare voids pertinent to 
these efforts. 

An interesting comparison can be made between our analysis 
of hydrogen reionization and similar effects during and after he- 
lium reionization (Dixon & Furlanetto, in preparation). In the latter 
case, the PDF of J is much broader (spanning over an order of 
magnitude), and the fluctuations have a much larger effect on the 
Hell Lya forest, introducing a factor ~ 2 bias in estimates of (r) 
when the background is assumed uniform. This has important im- 
plications for the Hell Lya forest during reionization: unlike the 
case we have studied here, the large fluctuations in the ionizing 
background lead to substantial opacity gaps even during the middle 
phases of helium reionization. 
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